Gravitational waves in dynamical spacetimes with matter content in the Fully Constrained 

Formulation 
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The Fully Constrained Formulation (FCF) of General Relativity (GR) is a framework introduced as an alter- 
native to the hyperbolic formulations traditionally used in numerical relativity. The FCF equations form a hybrid 
elliptic-hyperbolic system of equations including explicitly the constraints. We present an implicit-explicit nu- 
merical algorithm to solve the hyperbolic part, whereas the elliptic sector shares the form and properties with 
the well known Conformally Flat Condition (CFC) approximation. We show the stability and convergence pro- 
perties of the numerical scheme with numerical simulations of vacuum solutions. We have performed the first 
numerical evolutions of the coupled system of hydrodynamics and Einstein equations within FCF. As a proof 
of principle of the viability of the formalism, we present 2D axisymmetric simulations of an oscillating neutron 
star. In order to simplify the analysis we have neglected the back-reaction of the gravitational waves (GWs) 
into the dynamics, which is small (< 2%) for the system considered in this work. We use spherical coordinates 
grids which are well adapted for simulations of stars and allow for extended grids that marginally reach the 
wave zone. We have extracted the GW signature and compared to the Newtonian quadrupole and hexadecapole 
formulae. Both extraction methods show agreement within the numerical errors and the approximations used 
(~ 30%). 

PACS numbers: 04.25.D-, 04.30.Db, 04.40.Dg 



I. INTRODUCTION 

Numerical relativity is a rather young branch of physics 
devoted to the numerical solution of Einstein equations for 
complex problems, mostly in theoretical astrophysics, which 
involve the evolution of spacetime and eventually the matter 
content of a system. It was born thanks to the theoretical ad- 
vances which led to the 3+1 split of the Einstein equations 
IH 01 popularized due to the work of The 3+1 split de- 
fines a foliation of spacetime which allows to solve the equa- 
tions as an initial boundary value problem for a given spatial 
hypersurface which is then evolved in time. 

Soon after that theoretical breakthrough, the first hydro- 
dynamic calculations of the general-relativistic collapse of a 
star in spherical symmetry using a Lagrangian code were per- 
formed Jj, |5[] . Multidimensional simulations had to wait for 
the development of an Eulerian formulation [6] which could 
overcome the problems of non-spherical Lagrangian codes. 
Those advances led to a tremendously prolific era of general 
relativistic hydrodynamics with multiple applications to the 
formation of black holes, accretion onto compact objects, bi- 
nary neutron star mergers and core collapse supernovae (see 
[7] for a recent review on the topic). 

In parallel, a considerable effort was made to solve Einstein 
equations in vacuum. Only recently it has been possible to 
gain the sufficient understanding of the stability properties of 
numerical solutions of Einstein equations to overcome the nu- 
merical challenge of simulating the merger of two black holes 
and estimate its gravitational wave (GW) signal |8y. Despite 
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the fact that some of the first successful simulations 
used the Generalized Harmonic formulation (GH) 
of Einstein equations (see more about GH below), most of 
the groups in the numerical relativity community make use 
of a different numerical recipe, resulting from the combina- 
tion of different ingredients: i) the so-called BSSN formu- 
lation (JH [3; ii) the appropriate choice of gauge, with a 
slicing of the 1 +log family and some variant of the hyperbolic 
Gamma-driver condition for the spatial gauge Il5l4l7ll : iii) the 
use of high order spatial methods (at least fourth-order); and 
iv) the use of high resolution due to the increasing compu- 
tational power and development of adaptive mesh refinement 
(AMR) techniques. Using this recipe a number of groups 11181 - 
12311 have provided waveforms of one of the most powerful 
sources of gravitational radiation in the Universe, the binary 
black hole merger. These events are one of the main candi- 
dates for the first direct detection of GWs in the ground-based 
GW observatories (LIGO, Virgo, GEO600, TAMA300). 

According to the success of present formulations of the Ein- 
stein equations in different multidimensional scenarios, one 
could think that all the goals of numerical relativity have been 
achieved. Indeed, according to [24], the Holy Grail of numeri- 
cal relativity is a numerical code to solve Einstein equations, 
that simultaneously avoids singularities, handles black holes, 
maintains high accuracy, and runs forever. Current codes sat- 
isfy all the above requirements (with the exception, perhaps, 
of the last one), thanks to both the use of accurate numeri- 
cal techniques and stable formulations of Einstein equations. 
There are, however, a number of problems that could arise 
when using such an homogeneous set of numerical recipes to 
solve a complex multidimensional numerical problem. 

The first set of problems could be related to the homoge- 
neity itself. All working formulations of Einstein equations 
consist in the solution of hyperbolic equations, and in most 
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of the cases make use of the 3+1 decomposition of Einstein 
equations in a BSSN-like form, and use high order finite diffe- 
rences techniques to solve them numerically. This might make 
it difficult to detect fundamental theoretical problems in the set 
of equations to be solved. 

Some work has appeared facing these questions. The orig- 
inal (as far as solving the binary black hole problem) article 
of (8|] and some later ones flSJ UOll make use of GH formula- 
tion. It is a 4D covariant formulation which differs substan- 
tially from the 3+1 BSSN approach and allows for genuine 
comparisons between waveforms using different formulations 
l25l l26ll . Within the 3+1 decomposition of Einstein equa- 
tions, some formulations have appeared lfl5l l27ll in order to 
try to improve some of the weaknesses of the BSSN formu- 
lation, such as a better preservation of the constraints. An al- 
ternative is the characteristic approach (see B28I1 for a review). 
The Cauchy-characteristic matching and extraction, based on 
a (2+1)+ 1 formulation, has been successfully used to accu- 
rately extract gravitational waves matching its evolution to in- 
terior Cauchy data I29W31I1 . However, the success in using 
this approach to simulate the whole spacetime is very limited 
Jl^O due to the formation of coordinate caustics. Regard- 
ing the numerical methods, finite differences techniques have 
been used in all works, except for j9[] in which pseudo-spectral 
methods were used, although no substantial differences have 
been found in comparison with finite differences codes l26ll . 
Nevertheless, an alternative to pure hyperbolic formulations 
of the Einstein equations, as the work that we present here, is 
interesting and desirable. 

The second class of problems is related to the increasing 
level of complexity in the astrophysical scenarios that the 
community wants to simulate. After the binary black hole 
problem has been solved, the numerical relativity community 
is now concentrating more on the problem of solving non- 
vacuum spacetimes. The collapse of stellar cores and the 
merger of neutron stars represent challenges by itself, beyond 
the numerical evolution of Einstein equations: realistic micro- 
physics, accurate multidimensional neutrino transport, mag- 
netic fields, small scale instabilities (e.g., magnetorotational 
instability) and turbulence, non-ideal effects, elastic proper- 
ties of the crust of neutron stars, superfluidity and supercon- 
ductivity in cold neutron star interiors. Although full General 
Relativity (GR) is unavoidable in the case of the presence of 
black holes, approximations to GR in scenarios in which only 
neutron stars are present have a chance to simplify the nu- 
merical simulations to be able to understand the full physical 
complexity of those systems without the burden of solving the 
full GR equations. 

The typical neutron star has a mass of about M ~ 1 AM Q 
and a radius ofR ~ 10-15 km, which results in a compactness 
of GM/Rc 2 ~ 0.2 < 1. This implies that a post-Newtonian ex- 
pansion of the gravitational field of an isolated neutron star is 
possible and convergent, and that the expected error in the dy- 
namics of the system, if Newtonian gravity is used instead of 
full GR, is about 20%. However, it is well known that GR 
does not only produces quantitative effects in the dynamics of 
the system, but also qualitatively new effects like frame drag- 
ging or GWs. These two classes of effects appear at 1 and 2.5 



post-Newtonian level, respectively (see e.g. |H3l), which rep- 
resent changes in the dynamics of neutron stars of about 20% 
and < 2% at most. However, these effects can be important 
due to the non-linearity of the equations. As an example, in 
the case of neutron star binaries the energy loss due to GWs, 
despite their small nominal effect, make the orbit shrink until 
the neutron stars merge. In the case of the collapse of stellar 
cores, the stellar interior models used as progenitors can be 
treated safely in the Newtonian limit (GM/Rc 2 ~ 10~ 3 « 1), 
and only as nuclear density is reached GR effects appear. In 
that case, it would be desirable to have numerical tools which 
could allow us to evolve efficiently and smoothly the space- 
time, from the Newtonian regime of an stellar core to the 
mildly relativistic regime of a proto-neutron star or fully rela- 
tivistic regime of a black hole. 

An approximation to GR that could fill this gap between 
Newtonian gravity and GR is the Conformally Flat Condi- 
tion (CFC) approximation ( 135113611 ). The main features of the 
CFC approximation are: i) although it is not a post-Newtonian 
approximation, it behaves as a 1PN theory 113711 . and hence, it 
is possible to recover the Newtonian limit correctly in the case 
of weak gravity; ii) in spherical symmetry it coincides with 
GR, which makes it accurate for quasi-spherical objects like 
isolated neutron stars or for the collapse of stellar cores; iii) 
it only involves Poisson-like equations for the spacetime, and 
therefore the numerical methods and computational costs are 
closer to Newtonian simulations than to full GR simulations; 
and iv) it neglects GWs and the energy losses related to them. 
The numerical solution of elliptic equations is more involved 
than hyperbolic equations. However, the time-step in CFC is 
limited by the sound speed instead of the speed of light, as 
in the case of hyperbolic formulations of GR. That provides 
a considerably speed up in many scenarios that widely over- 
comes the extra cost of solving elliptic equations, making the 
numerical evolution considerably faster. The CFC approxi- 
mation was ori gina lly thought to deal with the neutron star 
binary case Il38l44lll . In this case, energy losses by GWs have 
to be included as an extra ingredient to allow for the neutron 
stars to merge. The major success, however, has been in the 
collapse of stellar cores ( Il42l444ll ). which lead to the computa- 
tion of GW emission using physically motivated microphysics 
j45tl . magnetic fields [46], and neutrino transport l47ll . The 
CFC approach has also been successfully used to simulate the 
phase-transition-induced collapse of rotating neutron stars to 
hybrid quark stars l48ll and the evolution of equilibrium mod- 
els of rotating neutron stars l49l l50ll . Direct comparisons of 
the CFC approach with full GR have shown that differences 
between both approaches, in the case of core collapse, are 
smaller than the numerical differences between the codes I5II - 
l53tl . This fact is understandable since the next post-Newtonian 
corrections to CFC were found to have an impact on the non- 
linear dynamics of less than 1% B54ll . In the case of neutron 
star mergers we are not aware of a direct comparison between 
CFC and full GR. 

A new formulation of Einstein equations which could ad- 
dress both classes of problems mentioned above, and share 
some properties with the CFC approximation, is the Fully 
Constrained Formulation (FCF) B55I1 . This formulation is 
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based on the 3+1 split of Einstein equations but different from 
all other formulations of Einstein equations that are purely 
hyperbolic, the FCF maximizes the number of elliptic equa- 
tions by solving the constraint equations at each time-step and 
choosing an appropriate gauge condition. As a consequence, 
the hyperbolic part of FCF only contains two degrees of free- 
dom, which correspond, far from the matter sources, to the 
GW content of the system. Therefore, FCF is fundamentally 
different from fully hyperbolic formulations of GR and can 
be used as another check of the consistency of the numerical 
solutions of Einstein equations. There are other formulations 
which incorporate the constraints into the evolution system 
(see 056[], section 5.2.2., for most relevant examples); how- 
ever, no numerical simulations have been performed yet with 
most of them (e.g., or simulations are restricted to axi- 

symmetric spacetimes (e.g., 115 811 ). Moreover, FCF is a natu- 
ral generalization of the CFC approximation; this fact makes 
possible a natural extension of all the numerical codes which 
use this approximation, in order to have a proper treatment of 
the gravitational radiation of the system without too much ef- 
fort. It also creates a bridge between weak gravity systems, 
which are well described within the CFC approximation, and 
the strong gravity limit. 

In practice, to extend a CFC code to FCF one has to add 
additional hyperbolic equations to the existing CFC elliptic 
equations (and also some extra sources in these elliptic equa- 
tions). This evolution system, written as a first-order one, is 
a hyperbolic system [59] and includes the whole hyperbolic 
sector of the metric of spacetime in this formulation. In partic- 
ular, the explicit values of the eigenvalues and eigenvectors of 
the hyperbolic metric system allow to guarantee the expected 
physical behavior on trapping horizons [60]. The remaining 
metric variables form the elliptic sector, which is similar to 
the group of elliptic equations in the CFC approximation, with 
extra sources. Recent works 16111 overcome some pathologi- 
cal problems related with non-local uniqueness in the elliptic 
equations in CFC, and also in FCF. The equations were rewrit- 
ten in such a way that these problems were solved, and the 
new scheme has been used successfully in some applications 




The analysis of the numerical evolution of the hyperbolic 
metric system is the main objective of this work, including 
aspects like numerical stability of the system in long-term 
simulations, evolution of equilibrium configurations, or the 
influence of the elliptic equations in the system. We per- 
form all the numerical simulations using the CoCoNuT code 
ll43l l44l l64ll . which was originally designed to evolve the 
hydrodynamics equations in the dynamical spacetime of the 
CFC approximation. The code uses spherical coordinates for 
the evolution of both matter and spacetime; this is very conve- 
nient in the present work, since it allows us to place the outer 
boundary sufficiently far from the star, in order to perform an 
accurate GW extraction 116511 . 

In this paper we present the first accurate extraction of the 
gravitational wave signature coming from the evolution of ro- 
tating oscillating neutron stars within FCF. As a first step to- 
wards a full evolution of the coupled system of elliptic and 
hyperbolic equations of the FCF, we have neglected the back- 



reaction of the GWs onto the dynamics of the system, which is 
a justified approximation in the case of isolated neutron stars 
and the collapse of stellar cores. 

The article is organized as follows. In Sec. II we review the 
FCF and detail the formulation used in the evolution of space- 
time. In Sec. Ill we describe the numerical methods used in 
the evolution of the different systems of equations (hydrody- 
namics, elliptic and hyperbolic metric equations). In Sec. IV 
we test our numerical implementation through the evolution 
of a vacuum spacetime with analytical solution. In Sec. V we 
perform simulations of equilibrium configurations of rotating 
neutron star and extract GWs from perturbed oscillating mod- 
els. Conclusions are drawn in Sec. V. Throughout the paper 
we use the signature (-, +, +, +) for the spacetime metric, and 
units in which c — G — M = 1 . Greek indices run from to 
3, whereas Latin ones from 1 to 3 only. 

II. SPACETIME EVOLUTION 

A. Fully Constrained Formalism 

Given an asymptotically flat spacetime (ALg,uv) we con- 
sider a 3 + 1 splitting by spacelike hypersurfaces denoting 
timelike unit normals to E, by rf. The spacetime on each 
spacelike hypersurface "L, is described by the pair {jij, K'i), 
where y^y - g MV + n^tiy is the Riemannian metric induced on 
Sf. We choose the convention K^ v - -jJLnJ^v for the extrin- 
sic curvature. With the lapse function N and the shift vector 
/?', the Lorentzian metric g^, can be expressed in coordinates 
(x?) as 

dxf dx v = -N 2 dt 2 + y/jidx' + ft dt)(dx j + /3 j dt). ( 1 ) 

As in 15511 we introduce a time independent flat metric fij, 
which satisfies -Ltfij = d,fij = and coincides with yy at 
spatial infinity. We define y := dety,y and / := detfij. We in- 
troduce the following conformal decomposition of the spatial 
metric y^y. 

y u = (//%■, if, = (r//) 1/12 - (2) 

The deviation of the conformal metric from the fiat fiducial 
one is denoted by h 1 ' , 

Hi : = fj - fJ. (3) 

Once the 3 + 1 conformal decomposition is performed, a 
choice of gauge is needed in order to properly reformulate 
Einstein equations. The prescriptions in 15511 are maximal 
slicing, 

K = 0, (4) 

and the so-called generalized Dirac gauge, 

Dif J = D t h ij = 0, (5) 

where K = y''Kij denotes the trace of the extrinsic curvature 
and Dk stands for the Levi-Civita connection associated with 
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the flat metric fy. More details can be found in [55]. Ein- 
stein equations then become a coupled elliptic -hyperbolic sys- 
tem: the elliptic sector acts on the variables if/, N, and ft 1 , and 
the hyperbolic sector acts on h'K More details of the analy- 
sis carried out for both elliptic and hyperbolic systems can be 
found in JH El. 

We introduce the conformal decomposition 

A' 7 := i[f m K ij , (6) 
and its decomposition in longitudinal and transverse-traceless 



parts 

A ij = (LX) ij + A^ T , (7) 

where 

(LX) IJ := D'X J + D J X' - -f' J D k X k (8) 

and DjA'^. - 0. These decompositions are motivated by 
the local uniqueness properties of elliptic equations shown 
in 06111 . We define w'i := T> k y 1 ' . The hyperbolic system for h! 1 
can be written as a first-order evolution system for the tensors 



~dt 



= 2Ni[f- b A ij + ffw[ ] - f k Dkp j - f'Dkp + -f j D k /3 k , 



(9) 



= £>k I ^-f'w 1 / + j3 k A ij ) - A kJ D k /3 l - A ik D k f3 j + ^A u D k j3 k + 2N^- 6 f k ,A ik A jl 



2 rll 



-SnNi/, 6 U 4 S ij -^- 



dw ' k 
dt 



N (if^R'J + Sy ik y jl D k if/Dnf/) + 4(A (y ik y jl D k if/DiN + y ik y jl D k NDnf/) 
~ [N[iff 2 R + 8y*'£>#£>^) + 8ifry kl D k if/DiN] y ij 
~\ {f kw k + f^^mW 1 ) - f k y ll D k D l (Nif) + ^y ij y kl D k D,W 2 ), 
D k \2Nif/- 6 A ij + p'w'/ - y u D,/3 j - y^Dip + \f j Dfi l 



where 



R = ^fDkhTDiymn ~ \fDkh mn D„y ml , 



(10) 
(ID 

(12) 
(13) 



S ij '■= Tftvyfy* is the stress tensor and S := 7 y Sy is its trace, T^ v being the energy-momentum tensor, measured by the observer 

of 4-velocity (Eulerian observer). Moreover, the system obeys the constraint of the Dirac gauge, w 1 . 1 = 0, and for the 
determinant of the conformal metric, we obtain y - f. 
The elliptic part of the FCF equations can be rewritten as 



u i . 7il7imA lm A^ iffR 



f'D k D,(M) = 



2nt//- 2 (E* +2S") + 



{lyjOjJM^ R\ 
8^ + ' 



y kl D k Di0+ \f k D k Di/3 l = \6nNif/- ( 'y i i(S*) j + A iJ Dj (2jV> -6 ) - 2Ni//~ 6 W kl A k 



(14) 
(15) 
(16) 



r 



The decomposition introduced in Eq. (|7]) leads to an extra 



where E := T MV n M n v and S , ■ :- -y^T flv n v are, respectively, the 
energy density and the momentum density measured by the elliptic equation for the vector X', 
observer of 4-velocity rf, E* := if/ b E, S* := i]/ b S, (S*)i := . . 1 

<//'.S\,. and 



DjSyr + -D'D k X k + y im [D k y ml - 



£>m?U 



(LX) k 



A u = 9^' { D 'f'j + ®j?n - D >yu) ■ (17) 



%7iy' ] (S*)j-y m \£> k y ml 



£> m yu 



TT' 



(18) 
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hydrodynamics 


elliptic metric 


hyperbolic 


Approach 


equations 


sector 


metric sector 


CFC 


[43] 


(23-CZ3) 


no 


Passive FCF 


[43] 




©-CD 


FCF 


[67] 




©-GD 


Teukolsky waves 


no (vacuum) 


fixed Minkowsky 


©-CD 


Equilibrium NS 


fixed, [72J 


fixed, [72J 


©-CD 


Oscillating NS 


[43] 




©-CD 



TABLE I. Guide to the approaches to GR used in this work: first 
three rows represent the approaches discussed in the theoretical part 
of the present work and the last three rows the approaches used in 
the numerical part. For each approach we provide the equations that 
we use for each of the three sectors (hydrodynamics, elliptic and 
hyperbolic metric sectors), a suitable reference where the equations 
can be found or a comment. 



and the evolution Eq. (TTOb can be viewed as an evolution equa- 
tion for the tensor A^ T . 

More details about the derivation of all the equations can 
be found in ll6lll . All these equations are to be solved coupled 
with the hydrodynamic equations for the evolution of matter 
which can be derived from the Bianchi identities and the con- 
tinuity equation, 



(20) 

(21) 
(22) 
(23) 



([T8T l. The resulting system of elliptic equations 

Ai/f = -mr'E* - fi,f ^ A '"' A '\ 

Si// 7 

A(Niff) =27rNi/r I (£* +2S*) 

+ Nil/- 7 -^-^ , 

v 8 

A0+ ^f ij DjD k j3 k = Dj(2Nilf- 6 A i3 ), 
AX 1 + ^DjD k X k = $nf ij (S*)j 



is identical to the CFC equations in the form described in 16111 . 

In the present work we solve the coupled evolution of the 
hyperbolic system for h'' given by Eqs. (l9]l-(fTTl). the elliptic 
approximated system for N, if/ and /? , given by Eqs. d20ll-(l2"2"l> 
and the hydrodynamics equations. We call the new system 
passive FCF, in the sense that we neglect the back-reaction of 
the GWs onto the dynamics of the system. Contrary to the 
CFC approximation, this approach does not neglect the GWs 
itself. Therefore, it allows one to compute the GW emission 
of the system directly from the spacetime evolution. Upper 
three rows of table U summarize the approximations used in 
the case of CFC, passive FCF and FCF. 



T% = 



4 = 0. 



(19) 



III. NUMERICAL METHODS 



Explicit expressions for the hydrodynamics equations for the 
case of a perfect fluid in the form that it is used in the present 
work can be found in 0641 • 



B. Passive FCF 



An interesting property of the fully constrained formalism 
is that if h'-* = 0, the resulting 3-metric y-,j is conformally 
flat. This condition corresponds to the well know conformally 
flat condition (CFC) approximation ll35l[36ll of Einstein equa- 
tions. The CFC approximation has been proved to provide ac- 
curate evolutions of spacetimes including single neutron stars 
and core collapse supernovae ll5ll - l53tl . In these scenarios the 
back-reaction of the GWs on the dynamics of the system is 
so small that h'i can be approximated to be zero. The main 
drawback of the CFC approximation is that the GW content 
is removed from the system, and the computation of the GW 
emission has to be performed approximately by means of the 
quadrupole formula. 

Since our aim is to deal with this kind of astrophysical sce- 
narios, neutron stars and core collapse supernovae, in which 
the GWs are not important for the dynamics, when solving the 
complete FCF system, the back-reaction of the h'i tensor onto 
the hydrodynamics and elliptic part of the metric equations 
can be neglected. Therefore, we impose h'-i = in Eqs. (fT4l- 



We perform all the simulations of this work using the nu- 
merical code CoCoNuT H |64j, HI] . We have extended this 
code, which solves the coupled evolution of the hydrodynam- 
ics equations with spacetime evolution in the CFC approxi- 
mation, to add the new degrees of freedom necessary for the 
FCF in the passive FCF approximation. In the following, we 
briefly describe the numerical methods used in the code to 
solve the hydrodynamics equations and the elliptic part of the 
FCF formalism. These methods and equations are identical 
to those described in ll6lll64ll . We also describe the numeri- 
cal techniques applied to solve the evolution of the h! 1 tensor, 
which is necessary to extend the CFC approximation to pas- 
sive FCF. In all cases we consider spherical coordinates and 
axisymmetry. In order to simplify the notation, we will re- 
fer to the three sets of variables as hydrodynamics variables, 
U := {D,Sj,r) (see definitions below), elliptic-spacetime va- 
riables or CFC variables, V := (N,tf/,f3',X'), and hyperbolic 
spacetime variables, W := (h lJ , A'i, w'J). 



A. Hydrodynamics equations 

The system of Eqs. ( fT9l can be cast into a system of con- 
servation laws [67] as 



<9,U + <3;F''(U,V) = <2(U,V). 



(24) 



U := (D,S j, t) is the conserved variables vector, D = —J^n^ 
and t = E — D. 
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Since we are neglecting the back-reaction of the GWs onto 
the dynamics of the fluid, there is no dependence on the 
hyperbolic-spacetime variables W in the previous set of equa- 
tions. We use Godunov-type schemes, which are suitable 
for solving equations written in conservative form. These 
schemes allow for a numerical evolution of the system with 
high accuracy in conservation of mass, momentum and en- 
ergy, and the correct behavior at discontinuities, e.g. shock 
waves at the surface of neutron stars (see e.g. 10]). We use 
the Marquina flux formula ll68ll combined with a second-order 
linear reconstruction with monotonized central slope lim- 
iter 16911 . The time update of the matter quantities relies on the 
method of lines in combination with a second-order accurate 
explicit Runge-Kutta scheme. The time step is restricted by 
the Courant-Friedrich-Lewi (CFL) condition [70|]. This com- 
bination provides second-order convergence in a number of 
tests including the evolution of oscillating neutron stars l64ll . 
which degrades to first-order in the presence of discontinu- 
ities. We use spherical coordinates (r,6,<p) in axisymmetry. 
The angular grid is equally spaced in 6, but the radial grid can 
be non-equidistant. 



B. Elliptic spacetime equations 

Once the values of the hydrodynamic variables, U, have 
been updated, the CFC metric, V, can be updated by solving 
the elliptic part of the spacetime evolution equations. It con- 
sists in a system of Poisson-like elliptic equations, Eqs. d20i>- 
d22l i. which can be written as 

AV = /(U,V). (25) 

This system of equations can be solved hierarchically, 
following the procedure described in l6lll . 

We compute the numerical solution using spectral methods. 
The sources of the equations are interpolated (parabolic inter- 
polation) from the finite difference grid to the spectral one, 
where the elliptic equations are solved using the LORENE li- 
brary for spectral methods 117 ill . The spectral solution of the 
equations is evaluated at the finite difference grid in order to 
update the metric fields N, 4> and /?, which are needed for the 
recovery of the primitive hydrodynamic variables (e.g., den- 
sity and velocity) and the evolution of W* (see next subsec- 
tion). The spectral grid consists of several radial domains in 
spherical coordinates. Further details can be found in |64]. 

Since the system of equations is elliptic, the Courant con- 
dition does not restrict the time step. Although the metric 
could be computed after every time step, in some scenarios 
the typical time scale of variation of the CFC variables is much 
longer than that of the hydrodynamic ones, and it is justified 
to compute V not at every hydrodynamics time step. In the 
simulations of neutron star oscillations presented in this paper 
we compute the CFC part of the metric every 10th hydrody- 
namical time steps and use a parabolic extrapolation between 
consecutive metric computations. This method has shown to 
provide sufficient accuracy for this scenario [43]. 



C. Hyperbolic spacetime equations 

Once we have updated the hydrodynamic variables, U, and 
the CFC variables, V, we solve the hyperbolic part of the 
spacetime evolution, Eqs. (l9l> — (fTTb. This part contains the 
gravitational wave information of the system. It consists of 
evolution equations for the variables W of the form 

d,W = s(W,V,U). (26) 

We solve the system following a two step approach. In 
the first step we update W' and A 1 ' to the next hypersurface, 
= Et»+A», denoted by an upperindex (n + 1), using only 
information of the previous hypersurface, denoted by («). 
It is therefore an explicit algorithm of the form 

d,h ij = S h (ti m , A ij(n \ wf n \ V w ), (27) 

d,A ij = S A (ti Kn \ A ij< "\ wf n \ U (n) , V (M) ), (28) 

which can be integrated using explicit Runge-Kutta schemes. 

In the second step we update wV using an implicit- 
explicit approach. We compute the sources using the values 
(h'M>,A'M>) and the updated values of (/^' ( " +1) , A* ,+1) ) com- 
puted in the first step. The sources of Eq. (TTTb can be splitted 
into two terms of the form 

d,w1 = S wl (ft", A' J , \) + S w2 (v^, V), (29) 

where 

S Kl = D k ^2Nifr~ 6 A ij - fD,/3 j - y'^fi + ^f J Dfij , 

S w2 = D k (p l w\ j ). (30) 

The first term, S w i, does not depend explicitly on the 
evolved variables w'/. The second term, S w2 , depends 
linearly on \vl and does not depend explicitly on the variables 
(h'j,A'i). This property allows us to design a numerical algo- 
rithm to evolve w'^ from E,» to , using the values of h'-* and 
A'-* at and all other variables at £j», i.e., 

d,Wl = S wX (h ij{n+x \A ij{n+l \ V (B) ) + S V)2 {wf n \ V (M) ). (31) 

This scheme provides a numerically stable evolution, due to 
the (partially) implicit dependence on S w \ and explicit on S w2 . 
We evolve the system with the same Runge-Kutta schemes as 
in the first step, but with the corresponding partially implicit 
evaluation of the S w \ source term (see Appendix lAl for more 
details). However, this evaluation reduces the theoretical or- 
der of the scheme, which is observed in numerical simulations 
in those scenarios where the terms S w i(h ,j( - n+1 \ A'j ( - n+1 \ V^) 
and S w i(ti i(n \A iiin \ V (n) ) differ significantly. In practice, the 
reduction of the order of the method can be small as long as 
the leading term of the sources for the evolution of wi is the 
one containing the A 1 ' tensor, as we have obtained in the evo- 
lution of Teukolsky waves with a method based on a fourth- 
order Runge-Kutta scheme (see Sec. II VI ). In other cases the 
order of convergence of the method can reduce up to second- 
order within the same scheme, as we have obtained when the 
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tensor h 1 ' reaches stationary values in the evolution of equili- 
brium configurations of rotating neutron stars (see Sec. lVBI ). 
In a general scenario, an implicit-explicit (IMEX) Runge- 
Kutta scheme lf73i p74h could be used to prevent the reduction 
of the order of convergence, although this is beyond the scope 
of this paper. The use of implicit terms for the second step 
of the time integration is crucial in order to provide stability. 
We have checked that when a purely explicit approach is used 



for wl, the numerical method becomes unstable 1 75] . The 
method becomes also unstable, when we compute w'J directly 
as spatial derivatives of h'K 

To solve Eqs. ( 1291 ) we use a fourth-order explicit TVD 
Runge-Kutta scheme B76I1 . together with the partially im- 
plicit treatment mentioned above. We use a fourth-order cell- 
centered Lagrange interpolation polynomials lf77ll to compute 
spatial derivatives, even for non-equidistan grids, and a forth- 
order Kreiss-Oliger dissipative term [78] to avoid the develop- 
ment of high frequency numerical noise. We impose an out- 
going radiation Sommerfeld 117911 condition to the linear part 
of the wave at the outer boundary, to prevent reflections from 
the boundary into the numerical domain. Unless other purely 
hyperbolic formulations of Einstein equations, the boundary 
conditions for the hyperbolic part of FCF do not influence 
the preservation of the constraints, since they are solved sepa- 
rately. 

The time step is determined by the Courant condition for 
the speed of light, c. This time step condition is more restric- 
tive than that of the hydrodynamics because the fluid eigen- 
values are bounded by c. The time step for W is chosen to be 
an integer fraction of the hydrodynamic time step, such that 
after each hydrodynamic time step U and W are synchronized. 



IV. TEUKOLSKY WAVES 

The first test is the evolution of Teukolsky waves [80] which 
are solution of the linearized Einstein equations in a vacuum 
spacetime. We choose as initial data a combination of ingoing 
and outgoing even parity axisymmetric Teukolsky waves with 
amplitude 10~ 5 . It provides regular initial data at r — which 
satisfies the Dirac gauge and is traceless (which is the linear 
approximation of unit determinant corresponding to orthonor- 
mal spherical coordinates for the conformal spatial metric y !; ). 
We keep the background flat, i.e., N - \f/ - I and f3' = 0. We 
assume symmetry with respect to the equatorial plane. The ra- 
dial interval [0, 10] and the angular one [Q,n/2] are discretized 
by n r and rig equally spaced grid points, respectively. Table U 
summarizes the approximations made in this test. 

We display in Fig.Q~|the radial profile of the component h n 
at the end of the simulation, t — 6, for three different numeri- 
cal resolutions n r x n$. Since the amplitude of the wave is 
sufficiently small to be considered a linear perturbation, we 
can compare the numerical solution with the analytical ex- 
pression for the Teukolsky wave at each time. The solution 
agrees in the propagation speed of the wave, its amplitude and 
the asymptotic decay with increasing radius with the analyt- 
ical solution. The Sommerfeld condition at the outer boun- 
dary produces ingoing reflections with at most an amplitude 
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FIG. 1. Radial profile of h" at / = 6 at the equator (upper panel) 
and at the pole (lower panel). Three different resolutions n T x n B are 
shown: 50x 10 (dotted lines), 100x20 (dot-dashed lines) and 200x40 
(solid lines). 



square of the outgoing wave, as it is prescribed by the imposed 
condition. In Fig. [2] we plot the absolute errors of the nu- 
merical solution with respect to the analytical values of all the 
non zero components of the tensor h 1 ' at (f, r, 6) = (6, 0, n/2) 
(Fig.Q]shows that the maximum of the absolute error appears 
at r — 0), for different resolutions. We obtain an order of con- 
vergence of 3.4, 3.6, 3.8, and 4.6 for the components h rr , h 09 , 
W and h' e , respectively, which is close to the fourth-order of 
the corresponding Runge-Kutta method. Note that, since the 
background has /3' = 0, the source term S W 2 — in Eq. (|29t ; 
in this case we observe no significant reduction of the conver- 
gence order. 



EVOLUTION OF EQUILIBRIUM ROTATING 
NEUTRON STARS 



To test the performance of the passive FCF formulation in 
an astrophysical scenario we perform simulations of the evo- 
lution of isolated neutron stars. In this case the gravitational 
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FIG. 2. Absolute errors of the numerical solutions with respect to the 
analytical values of all the non zero components of the tensor h' j at 
(t, r, 67) = (6, 0, n/2) in terms of the (scaled) cell size. Solid, dashed, 
dot-dashed and dot-dot-dashed lines fit the errors for the component 
h' r , h 00 , /z w and h , respectively. Dotted line is the reference of 
fourth-order of convergence. 
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FIG. 3. Scheme of the radial grid used in the code for rotating neu- 
tron star simulations. 
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field is sufficiently strong to need to go beyond the Newtonian 
limit and at the same time the GW back-reaction is sufficiently 
small for the passive FCF approximation to be valid. Table U 
summarizes the approximations made in the simulations of 
this section. 



A. Initial model and grid 

We construct the initial data using the numerical code rot- 
star jiirac of the lorene library [71], which computes axisym- 
metric and uniformly rotating neutron stars in equilibrium in 
the FCF formalism Il72ll . We use a poly tropic equation of state, 
P = Kp r , with r = 2 and K = 100 (in c = G = M = 1 
units), to construct a neutron star with a 550 Hz rotation fre- 
quency, an ADM mass Madm = 1.4874 M and a radius 
R = 15.18 km. The surface is located at a coordinate radius 
r, = 12.86 km at the equator. The wavelegth of GWs at the 
frequency of the fundamental f-mode (1.65 kHz as measured 
from our numerical simulations of Sec. IV CI ) is Ap = 181 km. 
Following IcUll the (local) wave zone for our neutron star 
model is located at r » A = A/(2n) = 28.8 km. 

Thanks to the use of spherical coordinates, we can adapt the 
radial grid to cover different domains of the space with the res- 
olution needed on each domain. In the case of a neutron star 
we have different resolution requirements inside the neutron 
star, where hydrodynamic variables have to be properly re- 
solved, and outside, where it is sufficient to resolve the wave- 
length of the outgoing GWs. 

For the finite difference grid, we consider three radial do- 
mains covering the computational domain (see Fig. [3): the 
matter domain (MD) contains the neutron star, extends from 
the center to a radius rMD-PD slightly larger than the stellar 
radius r t . This domain is covered by an equidistant radial 



TABLE II. Parameters used in the finite differences grid in simula- 
tions of neutron star evolution. 



grid. The propagation domain (PD), extends from ^d-pd to 
a radius tpd-dd >> In this region the radial grid spac- 
ing increases geometrically outwards, such that the GWs are 
well resolved. Near the outer edge of the domain, the GWs 
reach the wave zone, and hence it is an appropriate radius to 
perform the GW extraction. The damping domain (DD) ex- 
tends from ^d-dd to the outer boundary of the numerical grid 
'"out- We locate the outer boundary such that an outgoing wave 
generated at the center at t = and traveling at the speed of 
light reaches the outer boundary at the end of the simulation. 
This configuration minimizes the effect of spurious numeri- 
cal artifacts at the outer boundary which where found in our 
preliminary work l65ll . In the damping domain the radial grid 
spacing increases geometrically outwards and the wavelength 
of the GWs is not well resolved. This produces an effective 
damping of the outgoing GWs which helps to reduce the effect 
of the outer boundary conditions. To construct the finite dif- 
ference grid we need to provide the values of rMD-PD, '"pd-dd, 
r out , the number of grid points inside of the MD, the grid spac- 
ing ArpD-DD at ?"pd_dd (which automatically fixes the number 
of points inside PD), and the cell radial spacing ratio between 
consecutive ones at DD (which fixes the number of points in 
this domain). We perform simulations with two resolutions, 
labeled regular and high, whose grid parameters are given in 
table mi 

For the spectral grid we use two resolutions labeled regu- 
lar and high. The regular resolution grid consists of 5 radial 
spectral domains covering the finite difference grid, 4 domains 
with 33 and 1 with 17 radial collocation points, and a com- 
pactified spectral domain from r oul to infinity with 17 colloca- 
tion points. The 9 direction is covered by 5 collocation points. 
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FIG. 4. Radial profiles of U' for numerical simulations of neutron 
stars. The upper panel shows four snapshots of the evolution of h" 
and the initial stationary solution (grey line). The lower panel shows 
the components h" (black lines), h m (blue lines), and h w (red lines) 
for the initial stationary configuration (dotted lines), and the final 
configuration at 10 ms (dashed lines for the regular and solid lines 
for the high resolution spectral grid). 



The high resolution grid consist of 13 spectral domains, 6 do- 
mains with 33 and 7 with 17 collocation points, and a com- 
pactified domain with 17 collocation points. The 9 direction 
is covered by 17 collocation points. 



B. Equilibrium neutron stars 

Before attempting to solve the full coupled evolution of 
spacetime and hydrodynamics in neutron stars, we perform 
simulations in a more simplified setting. We evolve the hy- 
perbolic sector of the passive FCF formalism in a fixed non- 
trivial background (non-vanishing N,/3', if/ and hydrodynamic 
variables) corresponding to the equilibrium configuration of a 
rotating neutron star. Therefore, we evolve stationary initial 
data for the variable vector W and keep the variable vector U 



and V fixed during the evolution of W. 

Our fiducial model has regular resolution concerning both 
the finite difference grid and the spectral one. The background 
model is computed with FCF gravity, but we recompute the 
elliptic part, V, at the beginning with the CFC approxima- 
tion, i.e., the background is evolved in the same way as it is 
in the simulations of the next section, where the vector V is 
also evolved in time. The consequences of this modification 
on the background metric are evaluated below. We evolve the 
system of equations for the vector W for 10 ms. Due to small 
numerical discrepancies between the initial data for W and the 
numerical stationary solution of the equations a perturbation 
in the vector W appears and propagates outwards. The upper 
panel of Fig. [4] shows four snapshots of the evolution of the 
h rr component compared to the stationary solution (grey line). 
Note that the perturbation reaches the outer boundary at about 
3 ms, well before the end of the simulation. This is due to 
the unphysical superluminical propagation of the wave in the 
damping domain, where its wavelength is unresolved. How- 
ever, due to the smallness of h' 1 at the outer boundary, about 
3 orders of magnitude smaller than at the center, spurious re- 
flections are not noticeable in the simulations. At the end of 
the simulation the outgoing wave leaves the numerical domain 
and an equilibrium configuration remains. We can compare 
this solution with the initial stationary data to look for nu- 
merical discrepancies. At the center the initial configuration 
is recovered within ~ 10% accuracy. For distances r/r, > 10 
(see upper panel of Fig. |4j there are larger deviations, and the 
components of h'> decay approximately as instead of r -3 
as in the initial stationary model. Since GWs are contained in 
the part of h 1 ' decaying as r _1 , the erroneous decay observed in 
the simulations at large distances will lead to a constant offset 
in the computed GW amplitude. This offset can be compara- 
ble to the amplitude of GWs produced by small perturbations 
(see ll65ll for more details). We discuss the possible causes for 
this behavior below. 

One important reason for the appearance of offsets is the 
accuracy of the numerical solution of the elliptic equations. 
We have performed simulations increasing the spectral grid 
resolution, but keeping the same regular resolution for the fi- 
nite difference grid. The spectral grid resolution affects the 
accuracy of the computation of the background V computed 
at the beginning of the simulation. The lower panel of Fig. [4] 
shows the final configuration of the diagonal components of 
h 1 ' for the regular and high resolution spectral grid, compared 
to the initial equilibrium configuration. The error at the center 
of the equilibrium configuration at the end of the simulation 
does not improve with respect to the regular resolution spec- 
tral grid. The erroneous decay of W* improves significantly 
for the h rr component, and we recover the correct r~ 3 decay 
in the whole propagation domain. However, the h" and h w 
components do not improve significantly. Therefore, there is a 
strong sensitivity of the h rr component on the spectral metric 
resolution, and hence on the accuracy of the computation of 
V, because the variables in V appear in the leading terms of 
the equations for W, Eqs. (f9b — (TTTT). However, spectral resolu- 
tion does not seem to cure the problems in h ee and h w which, 
as we show below, are related to other sources of inaccuracy 
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FIG. 5. Time evolution of ft'' at r/r, = 19 .44 and 6 = tt/2. Left (right) panel shows the diagonal (non-diagonal) components. We plot two 
different finite differences resolutions: regular (thin lines) and high (thick lines). 
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FIG. 6. Same as Fig. [5] but we compare simulations using the CFC approximation (thin lines) and the full FCF metric (thick line) for the 
background metric V. 



in the solution of V. Because of the better performance of the 
high resolution spectral metric we use this resolution for all 
simulations hereafter. 

To check the effect of the finite difference grid resolution 
we have performed simulations with the high resolution grid, 
keeping the high resolution spectral grid fixed. The finite dif- 
ference grid resolution affects the accuracy of the solution of 
the vector W, i.e. the evolution of h l] . In Fig. [5] we plot the 
time evolution of the components of the tensor h! 1 (left panel 
for the diagonal components and right panel for non-diagonal 
components) at r/r* = 19.44, for the two finite difference grid 
resolutions. This coordinate radius is close to the outer boun- 
dary of the propagation domain and the inner edge of the ex- 
traction region for GWs (see next section). About 1 ms after 
the beginning of the simulation, the outgoing wave reaches 
this radius visible in Fig. [5] as a sudden rise of all compo- 
nents of ft'-*. After the outgoing wave leaves the numerical 
domain the value of h'> does not settle down to the initial 



equilibrium value, but to an offset value, decaying as r . All 
components converge with finite difference grid resolution to 
an offset value. The offset of the component h rr cannot be ap- 
preciated in Fig. [5] since it is much smaller than in the other 
components. 

Apart from the numerical error in the computation of V due 
to a finite spectral grid resolution, the approximating of the 
vector V by the solution of the CFC equations instead of the 
full FCF elliptic equations might also introduce small errors 
in the vector V, which are sufficiently large to explain the ob- 
served offsets. Although we still cannot solve the FCF elliptic 
equations with CoCoNuT in a simulation with spacetime evo- 
lution, for the case of fixed vectors U and V considered in this 
section, we can take the full FCF solution for V computed by 
the initial data solver rotstar jdirac : We have performed simu- 
lations with the regular and high resolution finite difference 
grids and the FCF background metric vector V. The resolu- 
tion used for the spectral solver in the initial data generator, 
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rotstarjdirac, fixes the numerical accuracy of the vector U. 
Unfortunately, due to internal code limits of rotstar jdirac , the 
maximum resolution that we could achieve was 8 radial do- 
mains, 5 with 33 and 3 with 17 collocation points, and 17 
collocation points in the 9 direction. This resolution lays in 
between the regular and the high resolution spectral grids used 
for CFC metric computations. Fig. [6] compares the evolution 
of h'i close to the GW extraction radius, r/r* = 19.44, with 
the CFC and the FCF background. The offset in h m , h w and 
h re is reduced when the CFC approximation is removed and 
the background is computed with the full general relativistic 
FCF formulation. The offset in h rr increases, however, al- 
though this is expected since the resolution of the spectral grid 
is lower in the FCF case than in the CFC case. We observe no 
change in the offset of the h r,p component. 

We conclude that the main reason for the offset in H' at 
large distance from the source is the accuracy of the compu- 
tation of the vector V. Both the resolution of the spectral grid 
and the neglected terms in the elliptic part due to the passive 
FCF approximation are responsible for this loss of accuracy, 
affecting in each case different components of h'i . This de- 
fines the spectral grid resolution necessary for the simulations 
in the next section, but the use of the passive FCF approxima- 
tion still introduces an offset which cannot be removed. The 
resolution tests show that the finite difference grid is adequate 
to evolve W, and is not responsible for the offset. 

The order of convergence of the code is difficult to evaluate 
due to the fact that if we increase the finite difference grid res- 
olution the solution converges towards an offset solution and 
not to the equilibrium one, and at the same time this offset con- 
verges with the spectral grid resolution and is affected by the 
passive FCF approximation. However, we can check the time 
behavior of the stationary solutions after the outgoing wave is 
gone. In Fig. [5] the time evolution of all /i'^-components suf- 
fers a time-drift which is due to numerical diffusion. We fit 
the evolution of the /i fifl -component between t\ — 4 ms and 
f 2 = 10 ms to h ee (ti) + C ■ (t 2 - t\) p . The fitted value for 
the power p is 1.7 and 1.8 for the regular and high resolution 
simulations, respectively. If we consider that C ~ Ar p , we 
can also compute the power as p — log(C regn i dr /Ci 1 i g t 1 )/log2, 
which is 1.96. Therefore, the order of convergence is close 
to second-order, due to the mixture of implicit and explicit 
terms in the fourth-order Runge-Kutta scheme. If we apply 
this analysis to other components of h' J , we find a similar or- 
der of convergence. 



C. Perturbed equilibrium configuration of rotating neutron 

star 

The last test consists in the evolution of an oscillating neu- 
tron star with coupled hydrodynamics and spacetime evolu- 
tion, i.e. we evolve the coupled system for U, V and W in 
the passive FCF approximation. We use the regular resolution 
finite difference grid and the high resolution spectral grid. We 
initiate the oscillations by adding a small I = 2 velocity per- 
turbation (about 1 % of the speed of light) to the stationary 
initial data used in the previous subsection. Previous prelim- 



inary studies l65ll show that the gravitational radiation has to 
be extracted close to ^pd-dd but still inside the propagation 
domain, where the wavelength of the fundamental mode is re- 
solved by about 5 grid points in the regular resolution finite 
difference grid. 

For asymptotically flat spacetimes, as the one used in our si- 
mulations, and since the Dirac gauge tends to the TT gauge far 
from the strong field region 18211 . it is possible to compute the 
amplitude of the plus polarization of the GW, i.e. the real part 
of the ^4 Weyl scalar, detected by an observer at a distance R 
with an observation angle with respect to the rotation axis 
as 



h JR, 0, T) = lim 



W\R, 0, T) - h 6e (R, 0, T) 



(32) 



In axisymmetry, the cross polarization, /i x , vanishes. In 
general, one should compute outgoing null geodesies for each 
angle 9 and determine the observation angle © and the dis- 
tance R, and then use the numerical value of h!' at the extrac- 
tion radius to determine h+(R,&,T). Thanks to the equato- 
rial symmetry, the curves 9 = n/2 on L, are null geodesies, 
and will be observed at a distance R with inclination angle 
© = n/2. The distance to an observer located at a coordinate 
radius r can be easily computed as 



Jo 



y rr {r, 9 = n/2)dr ■■ 



f 

Jo 



\l/(r,6 = nl2fdr. (33) 



For our neutron star model, the spacetime surrounding it is 
not extremely curved, and radial null geodesies at other an- 
gles are approximately curves with constant 9, i.e. » 9. If 
we integrate the distance R as in Eq. (1331 at different angles 9, 
the difference between equator and polar axis is about 0.04%. 
Hence, for neutron stars we can safely compute h+ at any an- 
gle using the numerical value of h' J at the extraction radius r ex t 
as 

, fl h^(r ext ,9,t)-h ee (r ext ,9,t)R(r exl ) 

h+(R,e,T)*> — , (34) 

Z K 

where T = t - R/c is the retarded time. Note that Eq. 041 ) is 
approximate and valid only in the wave zone, i.e. for r ext >> 
X ¥ . 

In axisymmetry and with equatorial symmetry, the multi- 
polar decomposition of the radiation field is 1 8 1 ] : 

h + (R, 9, T) = I (AfJ (T)T E2 - 2 \9) + A E2 T E2A "(9) + ...) (35) 

We have extracted the amplitude of the I = 2 and 4 multi- 
poles, A E q and A E ^, respectively, from the numerically com- 
puted h+, at all radial points between 15r» and 19.44r*. For 
post-Newtonian sources, as in our case, we expect the ampli- 
tude of the multipoles to decrease with / [81]. We were unable 
to extract higher-order multipoles because of the smallness of 
the signal, which made the numerical extraction too noisy. 
Since the extraction radii are relatively close to the source, 
in the interval 6.7/1-8.71, there is a small radial dependence of 
the GW amplitude and phase on the extraction radius. This 
dependence corresponds to r~ p , p > 2, components in h'-i, 




T [ms] 

FIG. 7. Quadrupolar component of the GW extracted from simulations of an oscillating neutron star. The upper panel shows the extrapolated 
waveform at infinity (black dashed line) and the offset-corrected waveform extrapolated at infinity (black solid line). The offset is corrected 
using the average value of the waveform of a non-perturbed neutron star simulation (black dotted line). The lower panel shows the offset- 
corrected waveform at an extraction radius of r/r, = 19.44 (black dashed line) and extrapolated at infinity (black solid line), compared to the 
quadrupole formula (orange dashed line). 



where the leading term is r~ 2 and can be fitted to extract the 
waveform at infinity using a procedure similar to f83ll . To ex- 
trapolate both phase and amplitude, we proceed in two steps: 
i) First, we perform a least squares fit of the retarded time of 
each maxima in the waveform as a function of the extraction 
radius R to 



T max (R) = r max (oo) + 



R 



(36) 



We use the average value of C for all maxima, < C > / 'r. = 
0.2245 ms, to correct the phase of the waveform as 



rM = T(R) 



<c> 

R 



(37) 



ii) Once the phase is corrected, we fit the amplitude of the 
waveform at constant T(oo) as a function of R to 



A%(R, r(oo)) = Ag(oo, T-(oo)) + 

K 



(38) 



The resulting fitted value of A™ (°°, T(oo)) is the waveform 
extrapolated at infinity. To estimate the finite distance effects 
we present results of the waveform extracted at a finite dis- 
tance r ext = 19.44r» and extrapolated at infinity. 

An alternative approach to compute the GW amplitude is 
to use the post-Newtonian wave-generation formalism. This 
is possible if the sources allow for a post-Newtonian expan- 
sion, i.e. (v/c) 2 ~ M/r„ < 1. For slow-motion sources 18111 . 



for which r* « A, it is possible to write the amplitudes of 
the different multipoles as volume integrals over the matter 
sources. Truncating the integrals at the lower post-Newtonian 
level, i.e. with Newtonian sources, the quadrupolar compo- 
nent results in the well known quadrupole formula 18411 . In 
axisymmetry the quadrupole formula reads 



<E2 
1 20 



n d 2 { r 



(3z?-l)r 2 y/ydrd0d<p)\, (39) 



with z = cos 6. Using the continuity equation, 7^ = 0, one of 
the time derivatives can be removed analytically (cf. [34]), 



A E2 _ 



-2>v* e z Vl -z 2 ) ryftdrdOdip) \, (40) 



where v" = av' - j3'. The latter formula is more convenient 
from the numerical point of view, since only one time deriva- 
tive has to be evaluated numerically. We use it in this work. 
In equatorial symmetry, the hexadecapolar component is 



<E2 
MO 



VBTr d A 
l26d7 



7z 4 - 6z 2 + - r 4 ^drd0d(p) 



(41) 




T [ms] 

FIG. 8. GW extracted from simulations of an oscillating neutron star. Upper and lower panels show the quadrupolar and hexadecapolar 
component respectively. On each panel we compare simulations for regular (dashed lines) and high (solid lines) resolution finite differences 
grid. The offset-corrected and waveform extrapolated at infinity computed with the direct extraction method (black lines) is compared to the 
PN method (quadrupole and hexadecapole formulae, orange lines). 



In a similar way as for the quadrupole formula, it is possible 
to remove one time derivative using the continuity equation 
(cf. JMUl): 



2 V5^ d 3 



Am ~ 63 dt 3 



+v* fl ((3-7z 2 )zVl -z 2 ) 



r 3 yftdrdOdip)}. (42) 



In the upper panel in Fig. [7] we plot the time evolution of 
the quadrupolar component A£? computed with Eq. ( [34-b {di- 
rect extraction hereafter) extrapolated at infinity. The wave- 
form clearly shows a constant offset which is unphysical. The 
dotted line shows the GW corresponding to the model in the 
previous section with the same resolution, but fixed U and V, 
and no initial perturbation. In this case we do not observe any 
oscillation, since the star itself is not oscillating, but an offset 
appears of similar magnitude as in the oscillating case. We can 
use the value of the offset at each time from the non-perturbed 
simulation to remove the offset in the oscillating simulation 
by simple subtraction of both GW signals. 

In the lower panel of Fig.|7]we compare the offset corrected 
waveform to the result of the quadrupole formula. There is a 
remarkable good agreement in the frequency and phase. The 
Fourier transform of the waveform is shown in the upper panel 
of Fig. [9] Within the numerical frequency resolution of the 
Fourier transform of the waveform, about 0.02 kHz, we do 
not observe differences in the frequency between the direct 



extraction method and the quadrupole formula. That sets an 
upper limit of 1 % for the frequency difference in the funda- 
mental oscillation mode, fit = 1 .66 kHz. The phase difference 
between both GW extraction methods, estimated as the rela- 
tive difference in the retarded time at the maxima, is smaller 
than 1%. The corrected signal (solid line in Fig. O agrees 
with the quadrupole formula within 30%. Therefore, the only 
big discrepancy with the quadrupole formula is due to the er- 
ror committed in the computation of the vector V, using the 
passive FCF approximation and the spectral grid resolution. 
The quadrupole formula is an approximate formula, which 
is valid in the slow-motion post-Newtonian limit. The error 
in the formula should be of the order (v/c) 2 ~ M/r*, which 
for our system is ~ 17%. Therefore, the 30% discrepancy 
in the amplitude is compatible with the approximation error 
of the quadrupole formula. Note that the waveform extracted 
at finite distance suffers from additional errors. For example, 
at 19.44r* (black dashed line in Fig. |7]l the amplitude of the 
waveform differs about 30% from the extrapolated waveform 
at infinity and its phase about 5%. 

In the case of the hexadecapolar component A^ 2 (lower 
panel of Figs.[8]and|9]l, the direct extraction method shows an 
important contribution due to numerical noise. Note that the 
hexadecapolar component is about a 2% contribution to the 
total waveform h+. Therefore, this numerical noise appears 
because of errors in the evolution of W below 1%, which are 
expected in our simulations. In this case the hexadecapole 
formula provides a good estimate of the phase and the fre- 
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FIG. 9. Fourier transform of the GW amplitude of an oscillating 
neutron star. Upper and lower panels show the quadrupolar and hex- 
adecapolar component respectively. On each panel we compare si- 
mulations for regular (dashed lines) and high (solid lines) resolution 
finite differences grid. The Fourier transform of the offset-corrected 
waveform computed with the direct extraction method (black lines) 
is compared to the quadrupole formula (orange lines). 



quency (~ 5%), however, the error in the amplitude is about 
50%. Note that for the hexadecapolar component there are 
possible sources of error in both the direct extraction method, 
due to the smallness of the amplitude, and in the hexadeca- 
pole formula, due to the three numerical time derivatives that 
we have to perform. Therefore, it is difficult to disentangle 
which one is a better approximation to the waveform. Nev- 
ertheless, both methods provide a reasonable agreement, and 
we are confident that the amplitude computed with any of the 
two methods is a rather good order-of-magnitude estimate of 
the hexadecapolar component. 

To test the effect of the finite difference grid resolution, 
we compare simulations with regular and high resolution, 
and the same high spectral metric resolution. Fig. [8] shows 
the offset-corrected values of Afg and A£? for both resolu- 
tions. For the offset correction we use the corresponding regu- 
lar and high resolution simulations of the previous section. 
We compare with the post-Newtonian wave generation for- 
malism (quadrupole and hexadecapole formulae, PN method 
for short) for both resolutions. In both methods, direct ex- 
traction and PN method, we observe a damping of the wave- 
forms which is reduced when increasing the resolution. The 



PN method only uses information in the GW generation zone 
(r < r„), contrary to the direct extraction method, where the 
wave is propagated from the generation zone to the wave zone. 
However, in both cases the damping observed in the wave- 
forms is of similar magnitude. That means that the source of 
the damping must be caused by numerical inaccuracies in the 
region close to the star, but not in the propagation domain. 
In other words, we observe a numerical damping of the os- 
cillations of the star itself, but not of the waves during their 
propagation towards the GW extraction point. 



VI. SUMMARY AND CONCLUSIONS 

We have reviewed the fully constrained formalism, which 
is a natural generalization of the CFC approximation of GR, 
and we have expressed the system of FCF equations in a form 
suitable for numerical simulations. We have presented a nu- 
merical scheme to solve the FCF system using spectral me- 
thods for the elliptic part and finite difference schemes for the 
hyperbolic part. In the simulations presented here we have 
neglected the back-reaction of the GWs onto the dynamics, 
which we call passive FCF. This work focuses on the stability 
and convergence of the hyperbolic part of the FCF equations, 
since the stability issues of the elliptic part were considered 
by llrJlll . We have presented a fourth-order finite difference 
scheme to solve the system of hyperbolic equations that makes 
use of implicit relations, to provide the necessary stability of 
the algorithm. We have solved the equations in spherical co- 
ordinates and axisymmetry. 

We have performed convergence tests for the hyperbolic 
part, in which the gravitational radiation of the system is en- 
coded, using a simple vacuum test with known analytical so- 
lution: the Teukolsky waves. We have shown the stability and 
convergence of the numerical evolution which is consistent 
with the fourth-order convergence of the numerical scheme. 
We have performed the evolution of equilibrium neutron stars 
and checked that the numerical code is able to maintain such 
configurations in equilibrium keeping the hyperbolic part to 
an accuracy of second-order. We interpret the drop to second- 
order convergence in all our simulations with matter content 
as an inconsistency in our numerical scheme in this regime, 
due to the mixture of explicit and implicit terms in the Runge- 
Kutta scheme. In order to improve the order of convergence 
one should use IMEX methods to solve the system of equa- 
tions 11731 17411 . Although this approach is beyond the scope of 
this work, it may be considered in the future. 

Finally, we have performed simulations of the evolution 
of an oscillating equilibrium neutron star. We have ex- 
tracted the GW signature from the metric components in 
the wave zone. We have compared the results from the 
direct extraction method with calculations using the post- 
Newtonian wave generation formalism, namely the Newto- 
nian quadrupole and hexadecapole formulae. We found that 
the approximate quadrupole formula describes within ~ 30% 
error the quadrupolar component of the wave. This is consis- 
tent with its nominal post-Newtonian error (~ 20%). Similar 
good agreement of the quadrupole formula with BSSN simu- 
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lations was found by l3ll[87ll . We were able to extract the he- 
xadecapole component of the wave, although numerical noise 
is considerably larger than in the quadrupole component. The 
comparison with the Newtonian hexadecapole formula agrees 
in frequency, phase and order of magnitude, however the com- 
parison is limited by the numerical accuracy of both wave ex- 
traction methods. We found that the evolution of the hyper- 
bolic part of the metric is very sensitive to inaccuracies in the 
elliptic sector, resulting in offsets in the GW signature. The 
main sources for inaccuracies are the number of collocation 
points in the spectral solver for the elliptic part and the ab- 
sence of back-reaction terms. We were able to get convergent 
results increasing the spectral resolution, however, an offset 
still remained due to the passive FCF approximation. We con- 
clude that GWs back-reaction should be included in the future, 
as well as an improvement in the accuracy of the numerical 
solution of the elliptic equations, in order to remove these off- 
sets. 
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We note that all simulations that we have performed in this 
work make use of spherical orthonormal coordinate grids, for 
the whole computation domain, in axisymmetry. Our simula- 
tions are among the few multidimensional simulations of Ein- 
stein equations in spherical coordinates. In the context of 3+1 
formulations, some of the first simulations of black hole for- 
mation used spherical coordinates [88, 89], however the for- 
mulations used in those works were unstable leading to expo- 
nentially growing constraint violations. Although some work 
has been done to reformulate the BSSN equations in order 
to ease its evolution in spherical coordinates I90II92I1 . these 
reformulations have been only tested for ID numerical prob- 
lems. The use of multi-block methods in BSSN simulations 
11931 l94ll allows to make use of spherical wavezone grids to 
take advantage of topologically adapted grids, but they keep 
a cartesian grid in the central region. On the other hand, 
spherical coordinates are widely used in the null formula- 
tions (see JUt]), mostly in the context of Cauchy-characteristic 
matching and extraction (e.g. l29T - [3llo . although stand-alone 
characteristic formulations have still very few numerical ap- 
plications (e.g. fUEl). 



Appendix A: Runge-Kutta schemes and partially implicit 
evaluation of the S w2 term 

In this appendix we describe with more details the numeri- 
cal method used in the evolution of the variables W, which has 
a partially implicit evaluation of the S ,, 2 term in the evolution 
of the w l J tensor. They are based in the explicit Runge-Kutta 
schemes of second, third and fourth-order. We show the pro- 
cedure with the second and third-order ones. 

The optimal second and third-order Runge-Kutta 
schemes B95I1 of a general evolution equation in time t 
for the variable u of the form u, = L(u) are respectively 

u m = u" + At L(u n ), 

M » +1 = I (u n + u m + AtL(u m j) , (Al) 

and 



We think that the reason for the success of our simulations 
in spherical coordinates is twofold. First, we use an implicit- 
explicit algorithm to solve the system of hyperbolic equations, 
whereas we solve implicitly the terms in the equations leading 
to instabilities. Second, only two degrees of freedom of the 
system, the GWs, are evolved by means of hyperbolic equa- 
tions, while the rest are the result of the computation of elliptic 
equations. This main feature of FCF is possibly crucial to pro- 
vide the extra stability to the numerical algorithm. We are not 
sure, whether both requirements are indeed necessary to per- 
form stable simulations in spherical coordinates, or whether 
the implicit-explicit scheme gives rise to the stability of the 
numerical algorithm. It would be interesting to explore the 
behavior of purely hyperbolic formulations with our implicit- 
explicit algorithm in spherical coordinates. 



m (1) = u n + At L(u n ), 



M (2 > 
u n+l 



3 M » + I M (l) + I Afi(M (l) x 

4 4 4 

I M » + 2 M (2' + -AfL( M ( 2) ), 
3 3 3 



(A2) 



where u n = u(t n ) and u n+l is the numerical approximation for 
the value u(f + At). 

The corresponding methods used in order to evolve W 
based on the previous Runge-Kutta schemes are: 
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1 . Based on the second-order Runge-Kutta scheme: 2. Based on the third-order Runge-Kutta scheme: 

h am = h m + A t s h (ti Kn \A iKn) ,wf n \v ( - n) ), 
#jm = &m + At s A {h^ n \A^ n \ wf n \ u (B) , v w ), w f = w f n) + AtS wi (h m \A m \ v (M) ) 



k 



4 4 4 

h m+D = Km) + i A w) + Ia* s^U'* 1 ', wf \ v w ), ^ = 3 AW0 1^ + 1 ijm Am wUm w w 

ILL 4 4 4 ^ v ' ' * ' ' " 



,4*4 

I 

3 3 3 



1 

— / 

2 



+ ^Af 5 w2 (wf\ V w ). (A3) ^i/Cn+D = LaUW + ^,7(2) 



3 3 
+ ^AfS i (^ (2) ) i y(2) ,wf 2) ,U , " , .V""). 

1 fifiil 2 f.™ 2 



,,'7(2) ,,(») 



3 



+ -AfS w2 (w^,V w ). (A4) 
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